Fortran 90 implementation of the Hartree-Fock 
approach within the CNDO/2 and INDO 

models 



Sridhar Sahu 1 , Alok Shukla 2 

Department of Physics, Indian Institute of Technology, Bombay, Powai, Mumbai 

400076, INDIA 



Abstract 

Despite the tremendous advances made by the ab initio theory of electronic structure 
of atoms and molecules, its applications are still not possible for very large systems. 
Therefore, semi-empirical model Hamiltonians based on the zero-differential overlap 
(ZDO) approach such as the Pariser-Parr-Pople, CNDO, INDO, etc. provide attrac- 
tive, and computationally tractable, alternatives to the ab initio treatment of large 
systems. In this paper we describe a Fortran 90 computer program developed by 
us, that uses CNDO/2 and INDO methods to solve Hartree-Fock(HF) equation for 
molecular systems. The INDO method can be used for the molecules containing the 
first-row atoms, while the CNDO/2 method is applicable to those containing both 
the first-, and the second-row, atoms. We have paid particular attention to compu- 
tational efficiency while developing the code, and, therefore, it allows us to perform 
calculations on large molecules such as Ceo on small computers within a matter of 
seconds. Besides being able to compute the molecular orbitals and total energies, our 
code is also able to compute properties such as the electric dipole moment, Mulliken 
population analysis, and linear optical absorption spectrum of the system. We also 
demonstrate how the program can be used to compute the total energy per unit 
cell of a polymer. The applications presented in this paper include small organic 
and inorganic molecules, fullerene Ceo, and model polymeric systems, viz., chains 
containing alternating boron and nitrogen atoms (BN chain), and carbon atoms (C 
chain) . 
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Program Summary 

Title of program: cindo.x 
Catalogue Identifier: 
Program summary URL: 

Program obtainable from: CPC Program Library, Queen's University of Belfast, 
N. Ireland 

Distribution format: tar.gz 
Computers : PC's/Linux 

Linux Distribution: Code was developed and tested on various recent versions 
of Fedora including Fedora 9 (kernel version 2.6.25-14) 
Programming language used: Fortran 90 

Compilers used: Program has been tested with Intel Fortran Compiler (non- 
commercial version 10.1) and gfortran compiler (gcc version 4.3.0) with opti- 
mization option -O. 

Libraries needed: This program needs to link with LAPACK/BLAS libraries 
compiled with the same compiler as the program. For the Intel Fortran Com- 
piler we used the ACML library version 3.6.0, while for gfortran compiler we 
used the libraries supplied with the Fedora distribution. 

Number of bytes in distributed program, including test data, etc.: size of the 
tar file bytes 

Number of lines in distributed program, including test data, etc. : lines in the 
tar file 

Card punching code: ASCII 

Nature of physical problem: A good starting description of the electronic struc- 
ture of extended many-electron systems such as molecules, clusters, and poly- 
mers, can be obtained using the Hartree-Fock (HF) method. Solution of HF 
equations within a fully ab initio formalism for large systems, however, is com- 
putationally quite expensive. For such systems, semi-empirical methods such 
as CNDO and INDO proposed by Pople and collaborators are quite attractive. 
The present program can solve the HF equations for both open- and closed- 
shell systems containing first- and second-row atoms using either the INDO 
model or the CNDO model. 

Method of Solution: The single-particle HF orbitals are expressed as linear 
combinations of the Slater-type orbital (STO) basis set specified by Pople 
and coworkers. Then using the parameters prescribed for the CNDO/INDO 
methods, the HF integro-differential equations are transformed into a matrix 
eigenvalue problem. Thereby, its solutions are obtained in a self-consistent 
manner, using methods of computational linear algebra. 
Unusual features of the program: None 
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1 Introduction 



The linear combination of atomic orbitals (LCAO) method is one of the most 
common approaches for solving the Schrodinger equation for many-electron 
systems such as atoms, molecules, clusters, and solids. It consists of express- 
ing the single-particle orbitals of the electrons of the system as a linear com- 
bination of a known basis set, and then solving the mean-field equations such 
as the Hartree-Fock (HF) or the Kohn-Sham equations. This converts these 
integro-differential equations into a matrix eigenvalue problem, which is subse- 
quently solved using computational approaches from the linear algebra. If one 
intends to go beyond the mean-field to include the electron correlations effects, 
approaches such as the configuration-interaction (CI), coupled-cluster, or the 
Green's function based formalisms can be used. If N is the total number of 
basis functions used, the computational difficulty at the mean-field level scales 
roughly as iV 4 , which is the number of two-electron integrals needed to perform 
such calculations. For post mean-field correlated calculations, integrals need to 
be transformed from the basis-set atomic orbital (AO) representation to the 
molecular-orbital representation (MO), a process which scales as A^ 5 , while 
subsequent solution of the corresponding equations can be even more time 
consuming[lJ. Since N increases rapidly with the number of atoms (and hence 
electrons) in the system, therefore, for very large systems solution of even the 
mean-field equations can become computationally intractable. Therefore, it is 
always advisable to devise methods of electronic structure theory which aim 
at reducing the size of the basis set. 

Using the zero-differential overlap (ZDO) approximation developed by Parr [2], 
Pople and coworkers developed a series of semi-empirical methods for comput- 
ing the electronic structure of molecules such as the Pariser-Parr-Pople (PPP) 
model|3], the complete neglect of differential overlap (CNDO) method [4|5|6] . 
and the intermediate neglect of differential overlap (INDO) methodf?]. Of 
these, the PPP model is applicable only to 7r-conjugated systems, however, 
the CNDO and INDO models with suitable parametrization, are in principle, 
applicable to all molecular systems [8] ■ CNDO and INDO methods are a class 
of valence-electron models which utilize a minimal Slater-type orbital (STO) 
basis set for the representation of the valence orbitals[8]. Additionally, in the 
representation of the Hamiltonian, only one- and two-center integrals are re- 
tained, leading to a drastic reduction in the computational effort as compared 
to the ab initio calculations|8j. Therefore, the CNDO/INDO models share at- 
tractive feature of semi-empirical parametrization with the PPP model, and a 
spatial representation of the molecular orbital with the ab initio approaches [8]. 
And, unlike the PPP model, the CNDO/INDO methods can also be used for 
the geometry optimization of molecules [8j. Therefore, for large molecular sys- 
tems and clusters, for which the applications of fully ab initio approaches can 
be computationally intractable, the CNDO/INDO methods provide an attrac- 
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tive alternative for the theoretical description of their electronic structure. 

It is with possible applications to large molecules, clusters, and polymers in 
mind that we have developed the present computer program which imple- 
ments the CNDO-2/INDO methods as formulated originally by Pople and 
coworkers[8j. As per the original formulation by Pople and coworkers [8], the 
INDO method can be used for the molecules containing the first-row atoms, 
while the CNDO/2 method is applicable to those containing both the first- 
, and the second-row, atoms. The fact that the code has been written in a 
modern programming language, viz., Fortran 90, allows it to utilize dynamic 
memory allocation, thereby freeing it from various array limits, and resultant 
artificial restrictions on the size of the molecules. Thus our program can be 
used on a given computer until all its available memory is exhausted. The 
present computer program can perform restricted Hartree-Fock (RHF) cal- 
culations on closed-shell systems, and unrestricted-Hartree-Fock (UHF) cal- 
culations on open-shell systems. Additionally, it also allows one to compute 
properties such as the molecular dipole moment, Mulliken population analysis, 
and linear-optical absorption spectrum under the electric-dipole approxima- 
tion. Apart from describing the computer program, we also present several 
of its applications which include various small molecules, fullerene Ceo, and 
polymeric chains consisting of carbon atoms (C-chain), and alternating boron 
and nitrogen atoms (BN-chain). 

The remainder of the paper is organized as follows. In section [2] we briefly 
review the theory associated with the CNDO/INDO approaches. Next, in 
section [3] we discuss the general structure of our computer program, and also 
describe its constituent subroutines. In section [4] we briefly describe how to 
install the program on a given computer system, and to prepare the input 
files. Results of various example calculations using our program are presented 
and discussed in section [5l Finally, in section [6J we present our conclusions, 
as well as discuss possible future directions. 



2 Theory 



In this section we briefly review the theory associated with the CNDO/INDO 
methods. The detailed discussion on the topic can be found in the book by 
Pople and Beveridge[8j. Our discussion will be in the context of the UHF 
method, the corresponding RHF equations can be easily deduced from them. 
As per the assumptions of the UHF method, we assume that the i-th up- 
and down-spin orbitals are distinct, and are represented (say) as ip^ and 
tpi^ , respectively. We assume that these orbitals can be written as a linear 
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combination of a finite-basis set 

# } = E^Vv, (i) 

where M 's represent the basis functions in question, and the determination of 
the unknown coefficients amounts to the solution of the UHF equations. 
In the equation above, we have only stated the expressions for the up-spin 
orbitals, the case of the down-spin orbitals can be easily deduced. Assuming 
the Born-Oppenheimer Hamiltonian for the electrons of the system 

fc2 N e N n N e 7 „2 „2 

« = -^Ev?-EE^ + E^. (2) 

zm i=\ ,4=14=1 n Ai i>j'i3 



where the first term represents the kinetic energy of N e electrons of the system, 
the second term represents the interaction energy of those electrons with its N n 
nuclei, Za represents the nuclear charge of the A-th nucleus, R^i denotes the 
distance between that nucleus and the z-th electron, represents the inter- 
electronic distance, while m and e are electronic mass and charge, respectively. 
We further assume that the total number of up-/down-spin electrons is N a /Np, 
such that N a + Np = N e . Using the conjecture of Eq. [Din conjunction with 
the Hamiltonian above, one obtains the so-called Pople-Nesbet equations 

- eZS^dg = 0, (3) 

V 



where S^ v is the basis function overlap matrix, ef is the UHF eigenvalue of the 
i-th up-spin orbital, is the Fock matrix for the up-spin electrons defined 
by 

= V + El p xAH^) - PxAH^)}, (4) 

Act 



above h^ u represents the matrix elements of the one-electron part (kinetic en- 
ergy and the electron-nucleus interaction) of the Hamiltonian of Eq.[2], (/ii/\\a) 
represents two-electron Coulomb repulsion integral in the Mulliken notation 

(fjiv | Act) = J J dT^MVMVruMtfMZ), ( 5 ) 

and P^a and P\ a , are the up-spin and total density matrix elements, respec- 
tively, defined as 

^ = £c£>c£\ (6) 

i=l 
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and 




(7) 



Equations [4] through [3, define the UHF method without any approximations. 
Next we briefly describe the approximations involved in the CNDO/INDO 
methods, leading up to corresponding UHF equations|8j: 

(1) Only valence electrons are treated explicitly thus N e = N v , where N v 
represents the number of valence electrons in the system. 

(2) A STO basis set centered on the individual atoms of the system is em- 
ployed, with the basis functions of the form 



where n, /, m represent the principal, orbital, and magnetic quantum num- 
bers associated with the /i-th basis function, Yi m (9, <fi) is the real spherical 
harmonic, and the radial part of the basis function is given by 



where is the orbital exponent associated with the /z-th basis func- 
tion, and is atom specific. In the CNDO method a minimal basis set is 
employed for the first row atoms, while an augmented basis set consist- 
ing also of d-type functions is employed for the second-row atoms. The 
present implementation of the INDO method, which is restricted only the 
to the first-row atoms, uses a basis set identical to the CNDO method. 

(3) The effective one-electron matrix elements are called core integrals, 
and determined semi-empirically. Diagonal elements of the one-electron 
element are determined through various parameters such as electron 
affinity A^, and ionization potential J M of the atoms involved, while the 
off-diagonal elements (// ^ v) are determined by 

h^v — PabS[iv, (10) 

where A and B denote the atoms on which basis functions \i and v are 
located, (3° AB is a semiempirical parameter dependent on A and B, and 
is the overlap matrix element for basis functions \i and v. 

(4) For orbital orthonormalization purposes it is assumed that the basis set 
is orthonormal. 

(5) For the two-electron integrals {p\> | Act), following approximation is adopted 



^(r,M)=i%(0n«(M), 
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Kiir) = (2C,)" +1/2 (2n!)- 1 / : 



V"- 1 exp(-C u r), 



(9) 



{\w | Act) = S^Sxai^ I AA) 



(11) 



and this set of integrals is further reduced by assuming 



(Hfi | AA) = jab, 



(12) 
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where it is assumed that basis functions /x and A belong to atoms A 
and B, respectively. The value jab is computed using the s-type orbitals 
located on A and B . Thus, all the two-electron integrals, apart from 
one- and two-center integrals, are ignored. As compared to the CNDO 
method, the following one-center integrals of the type (/iz/|/iz/) are as- 
sumed nonzero in the INDO method. The values of these integrals are 
determined semiempirically through the Slater-Condon parameters. 

Once all the approximations listed above are implemented, the diagonal ele- 
ments of the Fock matrix for the CNDO/2 model are given by 

= + A,) + Y,(Pbb ~ Z B ) lAB - (f~ - 1/2) 7 aa, (13) 

/ B 

while the off-diagonal elements for both the CNDO-2 and the INDO are 

= 0ab3^ ~ P^Iab- (14) 

In the equations above, Z B represents the effective nuclear charge of atom B, 
and Pbb = Z^eB P^n * s the sum of those diagonal elements of the total density 
matrix which are centered on atom B. In the INDO method, however, one uses 
different expressions for the one-center diagonal and off-diagonal elements, 
given by 

K> = u w + £ [^aa(HAA) - p^x^x)] 

+ E ( P BB ~ Z B ) lAB , (15) 

and 

f« = (2P^ - p;j(HH - C(HH, (is) 

where p,,u G A. Above U^, and the one-center two electron integrals are 
obtained through J M , A^, and various Slater-Condon parameters [7]. The two- 
center off-diagonal elements of the Fock matrix for the INDO model are ob- 
tained through [T4l Once the Fock matrix is constructed, both for the CNDO/2 
and INDO models, one solves the eigenvalue problem for the up-spin Fock ma- 
trix 

v 

as well as the down-spin Fock matrix, using the iterative diagonalization tech- 
nique, to achieve selfconsistency. From the equations given above, it is easy 
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to deduce the expressions for F^ u , as well as the Fock matrix elements for the 
RHF case. 



3 Description of the Program 

Our computer code consists of the main program, and various subroutines and 
modules, all of which have been written in Fortran 90 language. Additionally, 
the program must link to the LAPACK/BLAS library, whose diagonalization 
routines are used by our program. In the following we briefly describe the main 
program, and each subroutine. 

3.1 Main program CINDO 

This is the main program of our package which reads input data such as atomic 
numbers of the atoms constituting the system, and their positions, from the 
input file. The program also calculates the number of valence electrons of the 
system under consideration, and the total number of basis functions needed. 
It dynamically allocates various arrays, and then calls other subroutines to 
accomplish the remainder of the calculations. Because of the dynamical array 
allocation, the user need not worry about various array sizes, as the program 
will automatically terminate when it exhausts all the available memory on the 
computer. 

3.2 Subroutine BASEGEN 

This subroutine generates various arrays containing information to the basis 
functions used in the calculations. This includes quantities such as principal 
quantum number (n), orbital angular momentum (/), magnetic quantum num- 
ber (m), orbital exponent associated with each STO type basis function 
defined in Eqs. and [9l Additionally, it also stores some semi-empirical data 
associated with the Hamiltonian such the /3° and various Slater-Condon pa- 
rameters. This routine is called from the main program CINDO. 

3.3 Subroutine FACTCAL 

The primary task of this subroutine is to generate the factorials of various 
integers. The factorials thus generated are stored in global arrays accessible 
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via the MODULE factorials. This subroutine is also called from the main 
program CINDO. 



3.4 Subroutine ASSOC_LEGNDRE 



This subroutine initializes the expansion coefficients which define associated 
Legendre polynomials of various degrees, needed to represent the angular part 
of the basis functions. The data is stored in global arrays through MODULE 
legendre . This subroutine is also called from the main program. 



3.5 function SS 



A very important quantity used in computing Hamiltonian matrix elements 
is the so-called reduced overlap integral between two basis functions (labeled 
a, and b) [8 J 



CO 1 _^ 

s(n a , l a , m, n b , l b , a, (3) = J f(^ + v) na {ii - v) nb exp[--(a + (3)^} 

i -i 

x exp[— -(a — 0)v]T(fjL, u)dfidu, (18) 



where 



la—m l^—m 

T(fi,u) = D(l a ,l b ,m) E E Q amu C lbmu (^ ~ l) m (l - v 2 ) m 

U V 

x (1 + - H*0* + vY m ~ u {ti - vY m - v . (19) 



Above (n a ,Z a ,m) and (n b ,l b rn) are the quantum numbers of two basis func- 
tions, Ci amu , D(l a ,l b ,m) etc. are coefficients associated with the angular part 
of the basis functions, and a = ( a R, and (3 = ( b R, where £ a , ( b are the ba- 
sis function exponents, and R is the distance between the atoms on which 
basis functions are located. If we define the so-called Yy\ coefficients defined 
through the relation 



la— m h—m 

E E c lamu c hmu {y?-iT{i-v 2 r 

U V 

x(l + H"(l - nu) v {n + v) n *~ m -\ii - „)«>— = E Y^iS, (20) 

i,j=0 
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we obtain the expression 

s{n a , l a , m, n b , l b , a, ff) = D(l a , l b , m) ^ Y ijX Ai [-(a + /3)]S i [-(a - /?)], (21) 

i,3 



where 

oo 

A k {p) = J x k exp(— px)dx, (22) 



and 



Bk(p) = J x k exp(— px)dx. (23) 



For the s functions ( l a — l b — m — 0), the reduced overlap integrals (cf. Eq. 
T8l) can be written as 



s(n , 0, 0, n b , 0, a, /3) = - J J (// + i/) n<l (/i - z/) nb exp[--(a + /?)//] 

i -l 

x exp[-^(a - fi)v\dp,dv. (24) 
If we define the so-called Z k \ coefficients through 

n a +n b 

(/i + z/TGu - = J2 Z kxl i k ^ +nb ~ k \ (25) 

fc=0 



we obtain 

s(n a ,0,0,n 6 ,0,a,/3) = - £ ^^[-(a + /3)]S n< , +n6 _ fc [-(Q t - /?)]. (26) 

fc=0 ~* - 



The task of this REAL(kind=8) function is to compute the reduced overlap 
integral as defined in Eqs. I2T1 and 1261 for a given pair of basis functions a and 
b. The input to this routine is all the basis function related information such as 
their quantum numbers, orbital exponents, and the distance between them. It 
performs these calculations by calling subroutines GETYCOEF, GETZCOEF, 
AINT, and BINT which described below. 
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3.6 Subroutine GETYCOEF 



The task of this subroutine is to compute these Y^a coefficients, for a given 
set of n a , rib, l a , h, and m as defined in Eq. [20l It achieves this goal by calling 
subroutines BINOMIAL and POL2MUL described below. 



3.7 Subroutine GETZCOEF 

The task of this subroutine is to compute the Zk\ coefficients, defined in Eq. 
l25l for a given pair of s-type basis functions. As in case of subroutine GETY- 
COEF, this routine also computes for these coefficients by calling routines 
BINOMIAL and POL2MUL. 



3.8 Subroutine BINOMIAL 

Using the binomial expansion, expression (ax m y n + bx p y q ) 1 can be expanded 
as 



where i, j, m, n, p, q, and I, are integers, x, and y are variables, and a, b, and 
Cjj's are constants. This subroutine computes these expansion coefficients c^'s 
for a given set of input values of a, 6, m, n, p, q, and I. It is called both from 
routines GETZCOEF and GETYCOEF. 



3.9 Subroutine POL2MUL 

This subroutine computes the coefficients of the product polynomial when two 
polynomial of the type J2i,j a ij x% y 3 ^re multiplied, i.e., 



The input to this routine are coefficients c% and b^i, while the output consists 
of Cij. The arrays meant for storing these coefficients are allocated dynamically. 



(aa; m i/™ + 6a;V)' = E c ^V, 



(27) 



(28) 



i,j k,l l,m 
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3.10 Subroutine AINT 



The value of the integral A k (p), defined in Eq. [22j can be shown to be 

k\ 

^ p»(k - p + 



fc+i li 

A k (p) = exp(-p) -7TT ' , TT7 - (29) 



Subroutine AINT uses this series to compute the value of A k (p), for a given 
value of k and p. 



3.11 Subroutine BINT 



The purpose of this subroutine is to compute integral B k (p), defined in Eq. 
23l for a given value of k and p. We use the following recursion relation to 
perform the task 



B k+l (p) = -A k+ M+ ( - ir p e * M 

Thus first a call is made to the routine AINT to compute all the A k (p) 's 
needed. Subsequently, the 5 fc (p)'s are generated using the recursion relation 
of Eq. | 



3.12 Subroutine REDO VINT 



This is a very important subroutine which evaluates overlap matrix elements 
S^v among the basis functions. It evaluates the reduced overlap integrals de- 
scribed above for each pair of basis functions, by calling the function SS, using 
a coordinate system in which the atoms corresponding to the basis function 
pair are located along the z-axis. Then by a call to the subroutine TRANS 
described below, it obtains the actual overlap integrals by transforming the 
reduced integrals from the special coordinate system, to the actual molecular 
coordinate system. The upper-triangle of the overlap matrix is stored in a 
one-dimensional array. 
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3.13 Subroutine TRANS xmgrace 

The formulas for reduced overlap integrals (Eqs. [18] and [24]) assume that the 
atoms on which the basis functions are centered are a distance R apart from 
each other along the z-axis. But in practice, the molecules may have any kind 
of orientation. Therefore, we need to transform the reduced overlap integrals 
computed using these formulas, to the real orientation of the molecule. This is 
achieved through a transformation matrix which depends upon angular mo- 
menta of the basis functions, as well as on the angles by which the z-axis 
should be rotated to align it with the real orientation of the atoms involv- 
ing the two basis functions. The task of this subroutine is to construct this 
transformation matrix, and then apply it to obtain the overlap integrals with 
respect to the molecular frame. 



3.14 Subroutine COUL_INT 

This subroutine calculates the Coulomb integrals jab (c/. Eq. [12]) needed for 
the construction of the Fock matrix. It can be shown that these integrals are 
proportional to the reduced overlap integrals discussed above. Therefore, this 
routine computes these integrals by calling the function SS, and stores the 
values (one per atom pair) in a two-dimensional array. 

3.15 Subroutine CORE_INT 

The aim of this subroutine is to compute the one-electron part of Fock ma- 
trix, referred to as core integrals, and discussed in section [21 The calculation 
of off-diagonal elements involves the use of the overlap matrix elements S^ u 
computed in the routine REDOVINT, discussed earlier. The semiempirical 
data needed for computing these matrix elements is also passed to this rou- 
tine through arguments. The upper-triangle of the one-electron part of the 
Fock matrix, along with the extended Hiickel Hamiltonian, are finally stored 
in separate one-dimension arrays, and constitute the output of this routine. 



3.16 Subroutine DIPINT 

The aim of this subroutine is to compute matrix elements of dipole operator 
over the basis set. This subroutine is called only if the linear-optical absorption, 
or permanent electric dipole calculations are desired. Standard formulas are 
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utilized to compute these matrix elements, and it is called from the main 
program CINDO. 

3.17 Subroutine SCF_RHF 

This subroutine solves the RHF equations for the system under considera- 
tion in a self-consistent manner, using the iterative diagonalization procedure. 
The arrays which are needed during the calculations are allocated before the 
calculations begins, and are deallocated upon completion. Before the first iter- 
ation, extended Htickel Hamiltonian is diagonalized to obtain a set of starting 
orbitals. Subsequently, the Fock matrix corresponding to those orbitals is con- 
structed and diagonalized. The process is repeated until the self-consistency 
is achieved. During the self-consistency iterations, subroutine DSPEVX from 
the LAPACK/BLAS library is used to obtain the occupied eigenvalues and 
eigenvectors. Obtaining only the occupied eigenpairs, as against the entire 
spectrum, leads to considerable savings of CPU time for large systems. How- 
ever, if the entire spectrum of eigenvalues and eigenvectors is needed, say, to 
perform optical absorption calculations, the Fock matrix is diagonalized using 
the routine DSPEV from the LAPACK/BLAS library, upon completion of 
the self-consistency iterations. Because the entire spectrum is obtained only 
after self-consistency has been achieved, it does not strain the computational 
resources too much. Apart from computing the RHF total energy, this sub- 
routine also calculates the total binding energy of the system, and, if needed, 
performs Mulliken population analysis as well. 

3. 1 8 Subroutine SCF_ UHF 

This subroutine is exactly the same in its logic and structure as the previously 
described SCF_RHF, except that the task of this routine is to solve the UHF 
equation for the system under consideration. Different Fock matrices for the 
up- and the down-spin are constructed and diagonalized in each iteration, un- 
til the self-consistency is achieved. Similar to the case of routine SCF_RHF, 
during the iterations only the occupied eigenvalues and eigenvectors are com- 
puted using the routine DSPEVX. The iterations are stopped once the total 
UHF energy of the system converges to within a user defined threshold. 

3.19 Subroutine PROPERTY 

This is a driver subroutine whose task is to read the converged SCF orbitals 
written onto the disk by the SCF routines, and then call other subroutines 
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meant for computing various properties of the system under investigation. It 
is called from the main program CINDO after the SCF calculations, provided 
the user has opted for one of the property calculations such as the permanent 
electric dipole moment of the molecule, or its optical absorption spectrum. 

3.20 Subroutine DIPIND 

This subroutine transforms the dipole matrix elements from the basis-set AO 
representation to the SCF MO representation, by means of a two-index trans- 
formation. Therefore, it uses the dipole matrix elements computed in DIPINT, 
and the SCF MOs as inputs. The transformed dipole matrix elements, which 
constitute the output of this routine, are used in the calculation of linear op- 
tical absorption spectrum of the molecule. This subroutine is called from the 
routine PROPERTY described above, if the user has opted for the optical 
absorption calculations. 

3. 21 Subroutine DIPMOM_ RHF 

This subroutine calculates the total net electric dipole moment component of 
the molecule under investigation for restricted Hartree-Fock case. It is called 
from the routine PROPERTY if the user has opted for the dipole moment cal- 
culation. It uses dipole matrix elements calculated in the subroutine DIPINT 
and the SCF MOs as input, and computes the permanent dipole moment of 
the system using a straightforward formula. 

3.22 Subroutine DIPMOM_ UHF 

The purpose and logic of this routine is the same as DIPMOM_RHF, except 
that it is used for the case when UHF calculations have been performed. This 
routine is also called from the subroutine PROPERTY. 

3.23 Subroutine SPECTRUM 

This is an important subroutine which calculates the linear optical absorption 
of the system, under electric-dipole approximation, assuming a Lorentzian 
line shape and a constant line width for all the levels. Thus, if this calculation 
is opted, in the input file the user needs to provide the line width, along 
with the range of frequencies over which the spectrum needs to be computed. 
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Additionally, the routine uses the dipole matrix elements over the MOs as 
computed in routine DIPIND, along with the RHF single-particle energies. 
The computed spectrum is written in an ASCII file 'spectrum.dat', which can 
be readily used for plotting using programs such as gnuplot[9] or xmgrace[10j. 
This subroutine is also called from the routine PROPERTY, if the user has 
opted for linear absorption spectrum calculations. 

3.24 Orbital and Charge Density Plotting Subroutines 

For the purpose of orbital visualization, our code offers several options to 
the user for plotting the MOs, and the corresponding charge density. It is 
accomplished through four subroutines, PLOT1DRHF, PLOT 2D RHF, 
PLOT1DUHF, and PLOT 2D UHF. 

The task of subroutine PL0T_1D_RHF is to compute and print out the 
numerical values of RHF MOs, or their charge densities, on a one-dimensional 
grid of points, whose direction and range is provided by the user. Output 
of this program is written in an ASCII file called 'orbplot.dat', and can be 
readily used for plotting by gnuplot(9j and xmgracefTO]. This routine is called 
from the routine PROPERTY, if the user has opted for it. To compute the 
numerical values of RHF MOs (or their charge densities) at different points 
in space, it uses the numerical values of basis functions computed at those 
points, by calling function BASFUNC. 

When a user is interested in obtaining a two-dimensional plot of the RHF or- 
bitals/charge densities in the Cartesian planes, subroutine PLOT_2D_RHF 
is called from the routine PROPERTY. The structure of this routine is also 
similar to that of PLOT1D RHF, except that for this case the orbital/density 
values are printed out with respect to the two cartesian coordinates of the 
plane. This routine also uses function BASFUNC to compute the numerical 
values of the orbitals/densities, and the output is also written in the file 'orb- 
plot.dat'. In order to facilitate contour plots of charge densities, the option of 
making logarithmic plots is also available. 

In case of open-shell UHF calculations, the corresponding plots of the up- and 
down-spin MOs are obtained through calls to subroutines PL0T_1D_UHF 
and PLOT_2D_UHF, and the output is again written in the file 'orbplot.dat'. 

3.25 Function BASFUNC 

It is a REAL(kind=8) function whose aim is to calculate the numerical value 
of a given basis function, at a particular point in space. Therefore, the input 



16 



to this function consists of the coordinates of the point in space with respect 
to the location of the basis function, (n, I, m) quantum numbers of the basis 
function, and its exponent (. This function is called from all the orbital/density 
plotting subroutines described above. 



4 Installation, input files, output files 

We believe that the installation and execution of the program, as well as prepa- 
ration of suitable input files is fairly straightforward. Therefore, we will not 
discuss these topics in detail here. Instead, we refer the reader to the README 
file for details related to the installation and execution of the program. Ad- 
ditionally, the file 'input_prep.pdf explains how to prepare a sample input 
file. Several sample input and output files corresponding to various example 
runs are also provided with the package. 



5 Results and Discussions 

In this section, we present and discuss the numerical applications of our re- 
sults. First we present the results on a number of molecules. Next, we apply 
our method to obtain the ground states of model polymeric systems C chain 
and BN chain. Finally, we present the results of our calculations of optical 
absorption in Buckminster fullerene Ceo- Wherever possible, we compare our 
results to those published by other authors. 

5.1 Molecular Systems 

In this section we present the results of our calculations on a variety of 
molecules, including fullerene Ceo- The aim of these calculations is to com- 
pare our results with those published by other authors[8], and also with the 
CNDO/INDO calculations performed using Gaussian 03[H], in order to check 
the correctness of our program. 

In table ffl we compare the total HF energies of several molecules computed 
by our program, to those computed using Gaussian 03[H]. We used the bond 
length of 0.74 A for the H 2 molecule as used also by SurjanfTJ]. For water 
molecule we used the geometry from Schaeffer et al. |13| . for formic acid from 
Schwartz et a/. [14], for borazane from Palke[l5j, and for fluoropropene from the 
work of Scarzafava et al.[lQ\. For Ceo, we considered a dimerized configuration 
with the Ih symmetry group and bond lengths 1.449 A and 1.397 A optimized 
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Table 1 

Comparison of the total Hartree-Fock energies (Ejjf) of several molecules obtained 
using both the CNDO/2 and the INDO methods using our program, to those com- 
puted using Gaussian 03[1 1] . All results are in atomic units. S/E inside the parenthe- 
ses imply staggered/eclipsed configurations. For the molecular geometries utilized in 
these calculations, refer to the text. 



Molecule 


-E/nr(This work) 


Eh f (Gaussian03) 




CNDO/2 


INDO 


CNDO/2 


INDO 


H 2 


-1.474625 


-1.474625 


-1.474625 


-1.474625 


H 2 


-19.868052 


-19.013606 


-19.868052 


-19.013606 


cis-HCOOH 


-45.305164 


-43.364618 


-45.305163 


-43.364618 


trans-RCOOR 


-45.301984 


-43.360996 


-45.301984 


-43.360996 


Borazane (E) 


-20.169898 


-19.567825 


-20.169897 


-19.567824 


Borazane (S) 


-20.172764 


-19.570726 


-20.172763 


-19.570730 


cis-fluoropropene (E) 


-52.763845 


-50.693901 


-52.763845 


-50.693901 


cis-fluoropropene (S) 


-52.762138 


-50.692238 


-52.762138 


-50.692238 


irans-fluoropropene (E) 


-52.761824 


-50.692176 


-52.761823 


-50.692175 


£rans-fluoropropene (S) 


-52.759748 


-50.690019 


-52.759748 


-50.690019 


C60 


-427.624631 


-412.293447 


-427.624631 


-412.293447 



by Shibuya and Yoshitani [T7]. The Cartesian coordinates for the carbon atoms 
of Ceo were generated using the computer program developed by Dharamvir 
and Jindal[18]. Thus, for all the cases illustrated in the table, the agreement 
on the total HF energies between our calculations and those obtained using 
Gaussian 03Q3] is excellent both for the CNDO/2 and the INDO methods. 



Next we turn our attention to the comparison of results for geometry opti- 
mization of a few closed- and open-shell molecules. In table [2] we compare the 
bond lengths optimized by our program to those reported by Pople et al.[8\ 
for several closed- and open-shell diatomic molecules. Again the agreement ob- 
tained between the two sets of calculations is excellent both for the CNDO/2 
and the INDO calculations. 

Finally, in table [3] we compare the molecular dipole moments and Mulliken 
populations of several heteronuclear diatomic molecules obtained by our code 
with those reported by Pople et al. [8J. Both for the CNDO/2 and the INDO 
calculations the agreement between our results and those of Pople et al. [8] 
is virtually exact. Thus, excellent agreement between our results with those 
of other authors, not just for HF total energy, but also for other properties, 
testifies to the essential correctness of our computer program. 
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Table 2 

Comparison of geometries optimized by our code to those reported by Pople et al. [8], 
for several small molecules. Calculations for all the molecules with doublet or triplet 
ground states were performed using the UHF method. 



Molecule 


State 


Equilibrium Length (A) 






This work 


Pople et ai[B] 






CNDO/2 


INDO 


CNDO/2 


INDO 


Li2 


ly+ 
^9 


2.179 


2.134 


2.179 


2.134 


B 2 


ll 9 


1.278 


1.278 


1.278 


1.278 


c 2 


ly+ 
^9 


1.146 


1.148 


1.146 


1.148 


N+ 


2y+ 
^9 


1.127 


1.129 


1.127 


1.129 


N 2 


^9 


1.140 


1.147 


1.140 


1.147 


0+ 


2 n 

1L 9 


1.095 


1.100 


1.095 


1.100 


o 2 


3y- 
^9 


1.132 


1.140 


1.132 


1.140 


NH 




1.061 


1.069 


1.061 


1.070 


OH 


2 n, 


1.026 


1.033 


1.026 


1.033 


BeH 


2y+ 


1.324 


1.324 


1.324 


1.323 


LiH 




1.573 


1.572 


1.573 


1.572 


BN 


3 S+ 


1.269 


1.269 


1.268 


1.269 


LiF 


*E+ 


2.161 


2.162 


2.161 


2.162 


HF 


*E+ 


1.000 


1.005 


1.000 


1.006 


BF 


X E+ 


1.404 


1.408 


1.404 


1.408 



5.2 Calculations on Lithium Clusters 



In this section we discuss the optimized geometries of small lithium clusters 
computed using our program. The number of computational studies of the 
electronic structure of small lithium clusters by other authors is too numerous 
to list here. We will mainly refer to the ah initio works of Ray et a/.|20|. 
Boustani et al. [21], Jones et al. [22] , and Wheeler et a!. [23] who studied clusters 
similar to the ones studied by us. Detailed computational studies of several 
large atomic clusters containing various atoms are in progress in our group, 
and will be published later. 
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Table 3 

Comparison of computed electric clipole moments and Mulliken populations of het- 
eronuclear diatomics molecules with the work of Pople et al.[8\ The first number in 
each category is the CNDO/2 result, while the second number represents the INDO 
result. 



A/Tnl pr*n 1 p 


Electric Dipole Moment (Debye) 


Mulliken Population 




This work 


Pople et al. [8] 


This work 


Pople et al. [8] 


TYTTJ 

INJti 


1.76/1.69 


1.76/1.68 


A AO /A AA 

0.08/0.09 


A AO /A AA 

0.08/0.09 




1.78/1.80 


1.78/1.79 


0.16/0.18 


0.17/0.18 


BeH 


0.67/0.65 


0.67/0.64 


0.14/0.14 


0.14/0.14 


LiH 


6.16/6.20 


6.16/6.20 


0.27/0.29 


0.27/0.29 


BN 


0.36/0.50 


0.36/0.50 


0.05/0.03 


0.05/0.03 


LiF 


7.91/7.87 


7.90/7.86 


0.56/0.58 


0.56/0.58 


HF 


1.86/1.99 


1.86/1.98 


0.23/0.27 


0.23/0.27 


BF 


1.31/0.87 


1.31/0.86 


0.15/0.15 


0.15/0.15 



5.2.1 Li 2 



Results of our calculations on the optimized geometry of lithium dimer for the 
closed-shell ground state were presented in Table El As is obvious from the 
table that our optimized bond lengths of 2.179 A(CNDO) and 2.134 A(INDO) 
for L12 are in perfect agreement with similar calculations performed by Pople 
et al.[8\. As far as the comparison with the experiments is concerned, both 
these results are significantly smaller than the measured value of 2.672 A|19|. 
Therefore, it will be of considerable interest whether, or not, the inclusion of 
electron correlation effects will improve the results. 

5.2.2 Li 3 

Geometrical configurations for a triatomic cluster can be broadly classified 
as: (a) linear, and (b) triangular. For homonuclear systems such as Li 3 , the 
possible triangular geometries can be further subclassified into: (i) equilat- 
eral triangle, (ii) isosceles triangle, and (iii) a triangle with unequal arms. Of 
course, the equilateral triangle geometry (D 3h ) is expected to undergo Jahn- 
Teller distortion to a lower symmetry configuration. Indeed, several density- 
functional theory (DFT) and ab initio correlated calculations have indicated 
that the isosceles triangle geometry (G^u) is the most stable configuration 
for Liq [20ll2ll22|23| . Our calculations were performed on the doublet ground 
state using the UHF method, and the results are summarized in table BJ We 
found that equilateral triangular configuration is energetically more favorable 
as compared to the Jahn- Teller distorted isosceles triangles, as well as equidis- 
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Structure 


CNDO Results 


INDO Results 




Bond Length (A) 


E HF (a.u.) 


Bond Length (A) 


E HF (a.u.) 


Linear 


1.461 


-1.8870412 


1.457 


-1.8819986 


D 3h 


1.932 


-2.0133753 


1.919 


-2.0066737 



Table 4 

Optimized CNDO and INDO geometries of Li3, and corresponding HF energies 
(Ehf). Calculations were performed on the doublet ground states using the UHF 
method. 



Structure 


CNDO Results 


INDO Results 




Bond Length (A) 


E HF (a.u.) 


Bond Length (A) 


E HF (a.u.) 


Linear 


1.186 


-2.9683366 


1.185 


-2.9590571 


Square 


1.617 


-3.3083610 


1.612 


-3.2976927 



Table 5 

Optimized geometries of Li4 clusters of various shapes obtained by CNDO and INDO 
methods, and corresponding HF energies [Ejjf)- Calculations were performed on the 
closed-shell ground state. 

tant linear configuration, both for CNDO and INDO models. Optimized INDO 
and CNDO geometries are in very good agreement with each other. The po- 
tential energy surface of the triangular configuration shows interesting features 
for both the models. We find that if the equal arms of the triangle are longer 
or shorter than the optimized bond lengths of the D 3h geometry presented 
in table [H the system does exhibit Jahn- Teller instability. For bond lengths 
longer than those of the D 3 h geometry, the distorted triangle has an angle 
less than 60° between the equal arms, while for bond lengths smaller than 
the optimized values, the corresponding angle is more than 60°. However, the 
global minimum was found for the D 3h geometry described in table [H As far 
as the ab initio correlated and the DFT calculations are concerned, most of 
them report the length of equal arms of the Civ geometry close to 2.8 A, and 
the angle between them in excess of 70° [20l[2T|22f23| . Therefore, it will be in- 
teresting whether the inclusion of electron correlation effects will improve the 
agreement between CNDO/INDO models and the ab initio results. 



5.2.3 Li 4 

For Li 4 clusters various geometries, ranging from linear to tetrahedral are pos- 
sible as investigated, e.g., by Ray et al. [20]. However, as reported by various 
authors, a rhombus structure is energetically most favorable. As a demonstra- 
tion of our code we compute the relative stability of three possible structures 
of this system namely, linear, square, and rhombus. 
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Results of our calculations are summarized in table [5l We found that, of 
the three possible structures considered, the square structure had the min- 
imum energy. The rhombus shaped structures have lower energy than the 
square structure for bond lengths in excess of 2.1 A. However, the energies 
of those structures was found to be higher than those of square structures 
reported in table [5l This result is similar to what we obtained for Li3 for 
which the equilateral triangular structure was found to be more stable than 
the isosceles triangular structure, both within the CNDO and INDO mod- 
els. Our result for Li 4 disagrees with those obtained by correlated ab initio, 
and DFT calculations|20 ; 21 22 23j which predict the lowest energy for the 
rhombus structure with its acute angle close to 50°. Additionally, all ab ini- 
tio calculations predict bond lengths significantly larger than obtained here. 
Therefore, it is of considerable interest to explore whether the inclusion of 
electron-correlation effects will bring our results in better agreement with the 
ab initio ones. 

5. 3 Ground state of polymers 

Our code can be used to study both the ground and excited state properties of 
oligomers of various polymers because they are nothing but finite molecules, 
ranging in size from small to large. However, in this section we demonstrate 
that our code can also be used to obtain the ground state energy/cell, in the 
bulk limit, for one-dimensional periodic systems such as polymers. Thus, it 
can be used, e.g., for the purpose of ground-state geometry optimization of 
polymers, which is what we demonstrate next. 

The energy per unit cell of a one-dimensional periodic system can be obtained 
using the formula 

E cell = lim AE(n) = lim (E(n + 1) - E(n)), (31) 

n— >oo n— »oo 

where E(n + 1)/E(n) represent the total energies of oligomers containing 
n + l/n unit cells. Thus, using this formula, for sufficiently large value of n, one 
can obtain the energy/cell of a polymer in the bulk limit, from oligomer based 
calculations. In what follows we show that value of E ceU converges quite rapidly 
with respect to n, even for polymers which have metallic ground states. For 
the purpose of illustration we consider two model polymers namely chains con- 
sisting of: (a) carbon atoms (henceforth C-chain), and (b) alternating boron 
and nitrogen atoms (henceforth BN-chain) . A C-chain consisting of uniformly 
spaced atoms will be metallic, which, as per Peierls theorem[24], is not al- 
lowed. Therefore, such a system is expected to dimerize leading to an insulat- 
ing ground state|24|. On the other hand Peierls theorem is not applicable to 
the BN-chain, which is a band insulator and isolelectronic with the C-chain for 
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Number of unit cells 

Figure 1. Energy per two-atom unit cell, of an undimerized C-chain, plotted as a 
function of the number of unit cells. The C-C bond length was taken to be 1.297A. 



a two-atom unit cell. In an earlier from our group, we had studied the ground 
state geometry of C- and BN-chains using a fully ab initio methodology both 
at the RHF and the correlated levels, and concluded that C-chain does indeed 
exhibit dimerization, while the BN chain prefers the uniform geometry [25]. We 
explore the ground state geometries of these two systems using our code. In 
order to take care of the dangling bonds, we terminate the ends of oligomers 
of uniform C-chain and the BN chain with two hydrogen atoms on the each 
end. The dimerized C-chain consisting of alternating single and triple bonds, 
on the other hand, is terminated by one hydrogen atom on the each end. 



First we examine the convergence of E ceU obtained using Eq. [31] with respect 
to the number of unit cells n. In Figs. [Hand [2] we plot AE(n) as a function of 
n, for uniform C- and BN-chains, respectively. In both the cases convergence 
with respect to n, for the two-atom unit cells, is quite rapid, and for n = 10 
the bulk limit has been achieved to reasonable accuracy. This is quite remark- 
able because the C-chain considered for this calculation is metallic because 
of uniformly placed atoms. The convergence is even more rapid for C-chains 
with dimerized geometry. 
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Table 6 

Comparison of our CNDO/2 and INDO geometries for uniform C-chain, and the BN 
chain, with our earlier ab initio RHF results[25]. 



Calculation 


Bond Length (A) 




C-Chain 


BN-Chain 


This Work (CNDO/2) 


1.297 


1.360 


This Work (INDO) 


1.300 


1.362 



□ 



T 




13 



Number of unit cells 

Figure 2. Energy per unit cell of boron- nitrogen chain, plotted as a function of the 
number of unit cells. The B-N bond length was taken to be 1.360A. 



Results of our calculations are summarized in tables Ehndd From these tables, 
the following trends are obvious: (a) CNDO/2 and INDO optimized geometries 
in all the cases are in good agreement with each other, and (b) optimized bond 
lengths obtained here are slightly larger than those obtained using the ab initio 
RHF method^- 

Additionally, the condensation energy per atom of the C-chain, defined as the 
difference in E ce u per atom of the optimized geometries in the uniform and 
dimerized configurations, are obtained to be 11.1 mHartrees/atom (CNDO/2), 
and 11.3 mHartrees/atom (INDO). These numbers are in reasonable agree- 
ment with the corresponding ab initio RHF value of 7.8 mHartrees/atom [25]. 
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Table 7 

Comparison of our CNDO/2 and INDO geometries obtained for the dimerized C- 



chain, with our earlier ab initio RHF results[25 



Calculation 


'"singie(A) 




This work (CNDO/2) 


1.390 


1.231 


This work (INDO) 


1.390 


1.227 


Abdurahman et al.f2^ 


1.360 


1.174 



5.4 Optical Absorption in Fullerene Cqq 



Since the discovery of the Ceo in 1985|26], the field of the electronic structure 
and optical properties of fullerenes has become one of the foremost research 
topics these days[27]. Therefore, as the last application of our code in this 
paper, we present the results of linear optical absorption calculations in Ceo, 
at the RHF level. For these calculations we utilized the same geometry of 
Shibuya and Yoshitani[T7], as was used for total energy calculations presented 
in section 15.11 In the CNDO/INDO models, with four basis functions per 
carbon atom, Ceo has 120 occupied and 120 unoccupied orbitals, with the 
ground state being a closed shell with the A g symmetry. The HOMO/LUMO 
exhibit nearly tt/tt* character, with a five-fold degenerate HOMO (h u ) and a 
three-fold degenerate LUMO (t lu ). Our HOMO-LUMO gap of 9.23 eV for the 
INDO calcualtions is in perfect agreement with that reported by Shibuya and 
YoshitanifTT]. 
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800 




E (a.u.) 

Figure 3. Linear absorption spectrum of Ceo, obtained from RHF calculations using 
the INDO model, plotted as a function of the photon energy (in atomic units). A 
line width of 0.02 a.u. was assumed. 



Next we present the linear optical absorption spectrum of Ceo computed by 
the INDO method under the electric-dipole approximation, in Fig. [3j Because 
the HOMO and LUMO orbitals have the same inversion symmetry (unger- 
ade), the HOMO— >LUMO transition is dipole forbidden leading to negligible 
absorption intensity in the low-energy regions, in complete agreement with 
the experiments[27j. We have intentionally plotted the spectrum with a rel- 
atively small line width to emphasize the fact that a number of transitions 
among various orbitals contribute to the linear absorption. Qualitative fea- 
tures of our computed spectrum, namely the occurrence of two broad bands 
with a number of subpeaks in the spectrum, are in good agreement with other 
theoretical calculations [28j. As far as the quantitative comparison with the 
experiments is concerned, it is a well-known fact that the HF method over- 
estimates the energy gaps significantly. Therefore, in future works we intend 
to carry out various levels of CI calculations to investigate the influence of 
electron-correlation effects on the linear absorption in Ceo- 



6 Conclusions and Future Directions 



In this paper we have described our Fortran 90 program which solves the 
HF equations for both the closed- and open-shell molecular systems using the 
semiempirical CNDO/2 and INDO models. To demonstrate the correctness of 
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our approach, we presented numerous test calculations on molecular systems 
for which CNDO/INDO results are known, and obtained essentially exact 
agreement. Additionally, we presented results on systems such as clusters, 
fullerene, and polymers to demonstrate the wide utility of our present program. 
The reason behind developing the present program is twofold: (a) to develop 
a code in a modern language such as Fortran 90 which can carry out dynamic 
array allocation and thus free the user from specifying and changing array 
sizes, and (b) to provide an open software which will be widely available to 
users which they can use and modify as per their needs. One could write 
programs to perform a change of basis on the Hamiltonian matrix elements 
from the basis set AO representation to the MO representation, and use the 
transformed Hamiltonian to perform correlated CI calculations. Additionally, 
one could also introduce an electric-field in the Hamiltonian to perform finite- 
field calculations to compute quantities such as static polarizabilities of various 
orders. 

The present version of our code is restricted to first-row atoms using the INDO 
method and up to the second-row elements using the CNDO/2 approach. It 
will be extremely desirable to extend these methods to elements further in 
the periodic table, preferably up to the transition metals. However, there are 
several versions of these models available for heavier elements such as the s-p-d 
INDO, ZINDO, and other methods[29]. Therefore, one could implement these 
methods in the present code which will allow the user to perform both INDO 
and CNDO/2 calculations on elments of second-row and beyond. 

Work along those directions is continuing in our group, and results will be 
published as and when they become available. 
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